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Abstract 

Signal averaging is the process that consists in computing a mean shape from a set of noisy signals. 
In the presence of geometric variability in time in the data, the usual Euclidean mean of the raw data 
yields a mean pattern that does not reflect the typical shape of the observed signals. In this setting, it is 
necessary to use alignment techniques for a precise synchronization of the signals, and then to average 
the aligned data to obtain a consistent mean shape. In this paper, we study the numerical performances 
of Frechet means of curves which are extensions of the usual Euclidean mean to spaces endowed with 
non-Euclidean metrics. This yields a new algorithm for signal averaging without a reference template. 
We apply this approach to the estimation of a mean heart cycle from ECG records. 
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1 Introduction 

In many applications (biology, medicine, road traffic data) one observes a set of J signals that have a 
similar shape. This may lead to the assumption that such observations are random elements which vary 
around the same but unknown mean shape. Signal averaging is then the process that consists in computing 
a mean curve which reflects the typical shape of the observed signals. This procedure generally amounts 
to find an appropriate combination of the data to compute an average shape with a better signal-to- noise 
ratio. In many situations, the observed signals exhibit not only a classical source of random variation in 
amplitude, but also a less standard geometric source of variability in time. Due to this source of geometric 
variability, the usual Euclidean mean of the raw data may yield a mean curve that does not reflect the 
typical shape of the signals, as illustrated by the following application. 



1.1 Signal averaging in ECG data analysis 

An important application of signal averaging is the estimation of a mean heart cycle from electrocardio- 
gram (ECG) records. An ECG signal corresponds to the recording of the heart electrical activity. It is a 
signal, recorded over time, that is composed of the succession of cycles of contraction and release of the 



heart muscle. Each recorded cycle is a curve composed of a characteristic P-wave, reflecting the atrial 
depolarization, that is followed by the so-called QRS complex which corresponds to the depolarization of 
the ventricles, and which ends with a T-wave reflecting the repolarization of the heart (see e.g. [GH06] 
for a precise description of an ECG recording). In this paper, we present results on two different data 
sets from the MIT-BIH database [GAG"*" 13] . In Figure [D^a), we display data from an ECG record of 
a healthy subject having no arrhythmia, and in Figure ^h) we display data from an ECG record of a 
subject showing evidence of signiflcant arrhythmia (note that, in all the flgures showing ECG data, units 
on the vertical axis are in millivolts). 

In the analysis of ECG data, it is generally assumed that the heart electrical activity repeats it- 
self. Therefore, during an ECG record, one classically considers that the heart cycle of interest remains 
approximately the same with every beat, and that it is embedded in a random white noise with zero ex- 
pectation that is uncorrelated with the mean shape to be estimated. After an appropriate segmentation 
of an ECG record, one observes a set of signals of the same length such that each of them contains a 
single QRS complex. The preliminary segmentation step is done by taking segments (of the same length) 
in the ECG record that are centered around the easily detectable maxima of the QRS complex of the 
beats. Identiflcation of these maxima can be done using statistical methods to identity local extrema in 
noisy signals Big06 , [GK95) or by applying appropriate digital filters to identify typical parts of the QRS 
complex [PT85] . In this paper, we used the approach in [BigOG] to identify local maxima and to segment 
the ECG record into signals of the same length containing a single QRS complex. After this preliminary 
segmentation, one thus observes signals with approximately the same shape. In the normal case, we 
obtained J = 93 signals displayed in Figure O while in the case of cardiac arrhythmia we obtained J = 12 
signals displayed in Figure [3l 

To estimate the typical shape of a heart cycle and to improve the signal-to-noise ratio, one may use 
the Euclidean mean of these signals. In the normal case, this leads to satisfactory results as shown in 
Figure [4] since the average signal displayed in Figure IHa) clearly reflects the typical shape of the heart 
beats displayed in Figure [31^b)-(g). However, in the case of cardiac arrhythmia, the electrical activity of 
the heart is more irregular. This can be seen in the shape of the heart beats displayed in Figure [3l[|b)-(g) 
which may vary significantly from one pulse to another. Due to noise and geometric variability in the 
measurements, a simple averaging step may cause a low-pass filtering effect that leads to a mean cycle 
that does not reflect the typical shape of the heart beats in the ECG record |RR86bl IRR86al ILJB03] . 
In the case of cardiac arrhythmia, the Euclidean mean causes a low-pass filtering effect in the shape of 
the QRS complex, as shown by the results displayed in Figure O Indeed, it can be seen in Figure [3] 
that, around the time point t ~ 0.45 (which corresponds to the beginning of the QRS complex) there is 
rapid transition between a flat region and the peak of the R wave. For the Euclidean mean displayed in 
Figure [5]^a), this transition is slower, and it can be seen that the average heart cycle in Figure [5][^a) is 
a convolution by a smooth kernel of the signals in Figure [3l[|b)-(g) around the time point t ~ 0.45. To 
obtain better results, it is necessary to use alignment techniques for a precise synchronization of the heart 
beats, and then to average the aligned data. 



1.2 Deformable models 

Assume for simplicity that the signals are observed on the time interval [0,1], and that they can be 
extended outside [0, 1] by periodicity. An alignment technique consists in flnding a time synchronization 
of a set of signals. To be more precise, deflne a deformation operator (j)Q parametrized by G as a 
smooth increasing function i;^^ : [0, 1] — t- M such that 

(j)-^{t) = (t)-e{t) for all t G [0,1]. 
In the paper, we shall consider the following families of deformation operators: 
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(a) 
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(b) 

Figure 1: (a) ECG recording during 60 seconds of a healthy subject having no arrhythmia, (b) EGG 
recording during 60 seconds of a subject showing evidence of significant arrhythmia. 



(a) 



(b) 



(c) 




(d) 



(e) 



(f) 



Figure 2: Normal case: six signals containing a single QRS complex out of J = 93 extracted from the 
ECG recording displayed in Figure [^a) . Units on the horizontal axis are arbitrary. 
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(d) 



(e) 



(f) 



Figure 3: Case of cardiac arrhythmia: six signals containing a single QRS complex out of J = 72 extracted 
from the ECG recording displayed in Figure [^b). Units on the horizontal axis are arbitrary. 




Figure 4: Normal case: (a) Euclidean mean of the J = 93 signals after segmentation of the ECG record. 
(b)-(g) Superposition of six signals containing a single QRS complex (dashed curves) with the Euclidean 
mean (solid curve). 




(e) (f) (g) 

Figure 5: Case of cardiac arrhythmia: (a) Euchdean mean of the J = 72 signals after segmentation of the 
ECG record, (b)-(g) Superposition of six signals containing a single QRS complex (dashed curves) with 
the Euclidean mean (solid curve). 

- Translation operators: 4>0{t) = t — 6 and 4>Q^{t) = (f>^g{t) = t + 6, for 6 gM. (p = 1) and all t G [0, 1]. 

- Non-rigid operators: (pg : [0,1] — t- [0,1] is diffeomorphism of [0,1] parametrized by some £ MP, i.e. 

a smooth increasing function with 0^(0) = and (j)g{i) = 1 (a general method for constructing 
non-rigid deformation operators is described in Section [2. 4p . 

Given /i,/2 : [0,1] — )• R and a family {(j)e)e£Rp of deformation operators, the problem of time syn- 
chronization of two signals is to find a gW such fi{4>e{t)) ~ f2{t) for all t G [0,1]. In ECG data 
analysis, the most widely used alignment technique is time synchronization using translation operators 
by temporal or multiscale cross-correlation, see [LJB031 ITIRll] and references therein. 

In the presence of variability in amplitude and shape in the data, it is reasonable to suppose that the 
signals at hand satisfy the following deformable (or perturbation) model: 

y/ = /(00;(t^)) + ^,(00;(t£)) + fxej, j = 1, . . . , J, and ^ = 1, . . . , n, (1.1) 

where Yj denotes the ^-th observation for the j-th signal and ti = ^ are equi-spaced design points in [0, 1]. 
The function / : [0,1] — t- M in model (jl.ip is the unknown mean shape of the signals. The are supposed 
to be i.i.d. normal variables with zero expectation and variance 1. They represent a standard white noise 
in the measurements. The Zj's are supposed to be independent of the e^'s and are i.i.d. realizations of 
a smooth random process Z : [0, 1] — t- M with zero expectation. Therefore, the Zj's represent smooth 
variations in amplitude of the signals around /. Finally, the O^s are assumed to be i.i.d. random variables 
in W with zero expectation. Hence, the random deformation operators (pg* represent geometric variability 
in time in the data. 

Note that the deformable model (jl.ip is well adapted to ECG data processing. The main type of 
perturbations related to the analysis of ECG data (see e.g. [LJB031 ITIRllj ) are the baseline wandering 
effect (a low- frequency signal), electromyographic (EMC) noise, and power-line interference which is an 
amplitude and frequency varying sinusoid. This source of additive noise is modeled in (II. ip by the zero- 
mean processes Zj. The physiological nature of the electrocardiographic signal also alters the recording 
from heartbeat to heartbeat in lag and duration. In the ECG signal, there are therefore variations in 
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time of the heart cycle from one beat to another. This makes the heartbeats look shorter or longer 
in duration. After the segmentation of an ECG record into signals containing a single QRS complex, 
this local variability in time is modeled in (jl.ip by the non-rigid deformation operators (pg*. Aligning 
heart beats using non-rigid deformation operators is an alternative to the cross-correlation method which 
is classical used in ECG data analysis. One of the goals of this paper is to show that using non-rigid 
operators to align signals obtained by the segmentation of an ECG record may significantly improve 
the shape of the mean heart cycle computed by averaging the data after a standard alignment using 
cross-validation. 

1.3 Prechet means of curves 

The problem of estimating / and the deformation parameters 6* in model (jl.ip has been studied in 
[BCll] using the following procedure. First, for each j = 1, . . . , J, smooth the data (Yj)^^-^ to construct 
an estimator fj : [0, 1] — M of / o (pg* . In this paper, this smoothing step is done either by low-pass 
Fourier filtering or by wavelet thresholding. In a second step, estimate simultaneously the deformation 
parameters 0*, j = 1, . . . , J by minimizing the following criterion 

(0i,...,0j)= argmin M(0i, . . . , (1.2) 
(6»i,...,0j)Geo 

where 

M(0i, ...,Oj) = jpJ^' (^fjicp^e, it)) -J |] {(P-e, (t))) dt, (1.3) 

and 

0o = {(0i,...,0j) G (wy, ei + ... + ej = o}. 

Finally, in a third step take 

as an estimator of the mean shape /. 

As explained in |BCllj . the estimator / can be interpreted as a smoothed Frechet mean of the observed 
signals. The Frechet mean |Fre48) is an extension of the usual Euclidean mean to spaces endowed with non- 
Euclidean metrics. We refer to [Afsllj and [Hucllj for recent overviews of this notion and its application 
to the analysis of random variables taking their values in non- linear manifolds. Note that this procedure 
does not require the use of a reference template to compute the estimators Oi, . . . ,6j of the deformation 
parameters. Instead, one can interpret j J2j'=i fj'i'P^g ., (*)) as a template that is automatically computed 

during the minimization of (jl.3p . and onto which one searches to align all the curves fj. Therefore, this 
approach avoids the problem of taking one of the observed signals as a reference template which can lead 
to poor alignment results for data with a low signal-to-noise ratio. It should be also remarked that the 
constrained set ©o (onto which the minimization of M{Oi, . . . , 6j) is done) reflects the assumption that 
the deformation parameters Oj in model (II. ip have zero expectation. The choice of this constraint is also 
related to identifiability issues in model (jl.ip . and we refer to |BC11| for a detailed discussion on that 
point. 

In [BCllj . the statistical properties of / and the 0j's have been studied in details in the asymptotic 
setting n — )• -|-oo and/or J — )• -|-oo. However, the numerical performances of Frechet means of curves 
have not been much studied in [BC11| . The goal of this paper is then two folds. First, we describe new 
algorithms to compute / and the ^j's in the case of translation and non-rigid operators. We also discuss 
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the choice of the smoothing estimators fj. Secondly, we show the benefits of the methodology proposed 
in [BCll] for signal averaging of ECG data. 



1.4 Previous work on signal averaging 

In statistics, the problem of estimating the mean shape of a set of curves that differ by a time transfor- 
mation is usually referred to as the curve registration problem. It has received a lot of attention in the 
statistical literature over the last two decades. Among the various methods that have been proposed, one 
can distinguish between landmark-based approaches see e.g. |GK92j . Big06|, and time synchronization 



to align a set of curves see e.g. |RL01j . |WG97j . |LM04j . |KG88| . [TIRllj . 

All these approaches are based either on the use of one the observed signals as a reference template 
(see e.g. [TIRllj for curve alignment with translation operators), or on the Procrustean mean which is 
a standard algorithm to estimate a mean shape without a reference template. The Procrustean mean 
is based on an alternative scheme between estimation of deformation operators and averaging of back- 
transformed curves given estimated values of the deformation parameters. To be more precise, assume 
that Yj : [0, 1] — t- R denotes a linear interpolation of the data (i^^)]L^. Using our notations, Procrustean 
mean consists of an initialization step fo = j Z]/=i which is the Euclidean mean of the raw data taken 
as a first reference template. Then, at iteration 1 < i < imaxj it computes for all 1 < j < J the estimators 
9j^i of the j'-th deformation parameter as 



= arg min / \Yj{<t>g^{t)) - fi-i{t)?dt 



and then takes fi{t) = i Ei=i ^jC'/'e . ,(i))) 

as a new reference template. This procedure is repeated until 
the estimated reference template does not change. Usually the algorithm converges in a few steps. 

However, in the case of a low signal-to-noise ratio, it is preferable to first smooth the data before 
alignment. In this paper, to minimize the criterion ()1.3p . we propose an alternative algorithm to Pro- 
crustean mean. We also show that the obtained smoothed Frechet mean ()1.4p yields better results. In 
the statistical literature, the criterion (11. 2p has been first proposed by [GLM07j in the case of translation 
operators and then further studied by [VimlOj and [BGlOj . Its generalization to other deformation oper- 
ators and the connection between minimizing ()1.3p and the computation of Frechet means of curves has 
been investigated in |BC11| . 



1.5 Organization of the paper 

In Section [2] we describe a gradient descent algorithm to compute the Oj^s. We detail this procedure in 
the case of translation and non-rigid operators. We also compare our approach with the Procrustean 
mean of curves. Numerical experiments on signal averaging for ECG data are presented in Section [3l We 
conclude the paper by a brief discussion on the numerical results and the benefits of our approach. 



2 Methodology for mean pattern estimation 

2.1 A gradient descent algorithm with an adaptive step 

Denote by 

VM{e) = {J^-Mio), . . . , JjM(0)) G {wy\ 

the gradient of the criterion ()1.3p ai = {Oi, . . . , 6j) G (M^)"^. A closed form expression of this gradient 
is given in Section [2.31 for translation operators, and in Section [2.41 for non-rigid operators. Computation 
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of the deformation parameters can then be done by minimizing the criterion M{0) using the following 
gradient descent algorithm with the constraint that % = 0: 

Initialization : let 0° = 0, 5o = fvATcF)!' ^^^^ = M(0°), and set i = 0. 



Step 2 : let = 6>* - 5^VM{e') and let 6'"'^"' = (^^i - i ^^^^ ~ej,...,~ej-^ Yi^i ^j) ■ 
Let C{i + 1) = M(6'"""'). 
While C(i + 1) > C{i) do 

and set C{i + l) = M(6>"""'). 
End while 
Take 0^+^ = 0"^^". 

Step 3 : if C{i) - C{i + 1) > p(C(l) - C{i + 1)) then set z = z + 1 and return to Step 2, else stop the 
iterations, and take 6 = 0*"''^. 

In the above algorithm, p > is a stopping parameter and k > 1 is a parameter to control the choice 
of the adaptive step 5i . 

2.2 Choice of the regularization parameter in the smoothing step 

For the smoothing step, we present numerical results for 

- low-pass Fourier filtering: for t £ [0, 1] 

|A:l<Aj 

with c"^^ = i Yll=i Yj/e~'^'^^''~ , and where \j S N is a regularization parameter (cut-off frequency). 
A possible data-based choice for Xj is to use generalized cross validation (GCV), see e.g. jCW79j . 

- wavelet smoothing by hard thresholding: for t G [0, 1] 



2^0 mi 2™ 

k=0 m=mo k=0 U ' I J 



™0 



where 4>mo,k{*) = 2—(l>{2"'H - k) and Vm,fc(i) = 2-^(2"* - k) are the usual scaling and wavelet 
basis functions at resolution levels < mo < m < mi and location k, k^f^m k respectively 
the empirical scaling and wavelet coefficients computed from the data (1^^)"^^ (see e.g. |ABS01| 
for further details on wavelet thresholding). The universal threshold aj^J2\og{n) depends on the 
estimation dj of the level of noise in the j-th signal. It is given by the median absolute deviation 
(MAD) of the empirical wavelet coefficients at the highest level of resolution mi (see e.g. [ABSOlj ). 
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2.3 The case of translation operators 

When (pg{t) = t — 9, it follows by simple computations that for j = 1, . . . , J 

-^MiO) = -| p,it + e,) Yl frit + 0,0 j + 1 [m + 0,)? - me,)?) . 

Note that in the case of low-pass Fourier filtering, fj{t + 1) = fj{t) for any t G M, and the gradient of M 
can be computed directly in the Fourier domain thanks to Parseval's relation 



de 



J 



\k\<\ 



J 

j'=l 



^(i')g2jfc7r6»j., 



where 3^(c) denotes the real part of a complex number c. 
2.4 The case of non-rigid operators 

To build a family {(j)9)e£M.p of parametric diffeomorphisms of [0, 1], we adapt to one-dimensional curves 
the approach proposed in [BGL09] to compute the mean pattern of a set of two-dimensional images. Let 
u : [0, 1] — 7- M be a smooth parametric vector field given by a linear combination of p basis functions 
{hk : [0, 1] M, k = 1,... ,p}, such that 



v{t)=Y,ekhkit) forte [0,1], 



k=l 

where 9 = {9i, . . . ,6p) G is a set of reals coefficients. The function v is thus parametrized by the set 
of coefficients 6, and we write v = vg to stress this dependency. In what follows, it will be assumed that 
the basis functions are continuously differentiable on [0, 1] and such that and ■§ihk vanish at t = and 
t = 1. For the /i^'s we took in our numerical experiments a set of p = 10 B-spline functions of degree 3 
using equally-spaced knots on [0, 1]. Then, let t £ [0, 1] and for u € [0, 1] consider the following ordinary 
differential equation (ODE) 

—i;{u,t)=ve{^l^{u,t)) (2.1) 



with initial condition ip{0,t) = t. Then, it can be shown (see e.g. [YoulO] ) that for any u £ [0,1] the 
solution of the above ODE is unique and such that 



t ^ il)0{u,t) = t + I ve{ipe{s,t))ds 
Jo 

is a diffeomorphism of [0, 1] i.e. a smooth increasing function with ipg(u,0) = and il^o{u, 1) = 0. Then, 
we denote by 



the solution at u = 1 of the ODE (12. ip . In this way, we finally obtain a diffeomorphism (f>0 that is 
parametrized by the set of coefficients 6 G M^, and that is such that (t>Q^{t) = (j)-0{t). Then, it can be 
also shown that the mapping 9 ^ (j)0 is differentiable, and thanks to Lemma 2.1 in |BMTY05] it follows 
that 

^e^t) = —J)0{t) \ \—J)0{u,t)\ hk{ilJe{u,t))du. 



dt 



dt 
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Now, one obtains that if (l)g is a diffeomorphism generated by the ODE (12. ip . then for 6j = {6f \ • • • , ^j^^) ^ 



with for g = 1, . . . ,p 

d 



(9) 



and 



5# ^ ^ 



2.5 Comparison between the smoothed Prechet mean and the Procrustes mean 

To illustrate the advantages of the smoothed Prechet mean over the Procrustes mean of the raw data 
(without any smoothing), let us consider a set of J = 15 signals generated from the following deformable 
model using translation operators: 

Yf = f{U - 0*) + Z.ite - e*) + ae], j = 1, . . . , J, and £ = 1, . . . , n, (2.2) 

where n = 128, the 0*'s are i.i.d. normal variables with zero mean and variance /x^ = 0.004, the e^'s are 
i.i.d. normal variables with zero mean and variance cr^, / is the signal displayed in Figure [B^a), and the 
Zj^s are i.i.d. copies of the Gaussian process Z defined by 

+ 00 

Z{t) =aQ + Y, k~^^HakCk{t) + &fcSfc(i)), t G [0, 1], 
k=l 

where qq and (a^, bk)k>i are i.i.d. normal variables with zero mean and variance cr^, Ck{t) = \/2 cos(27r/ct) 
and Sk{t) = \/2 sin(27rA:t). A sample of three signals out of J = 15 generated from model (12. 2p is displayed 
in Figure [6l^b)-(d). It can be seen that a has been chosen to have a relatively low signal-to-noise ratio. 
The estimation of / by smoothed Frechet mean using low-pass Fourier filtering with Xj chosen by GOV 
for each j = 1, . . . , J is displayed in Figure [6l|e) , while the estimation obtained by Procrustes mean of the 
data without smoothing is shown in Figure [IJf ) . Clearly, the result obtained by the smoothed Frechet 
mean is much more satisfactory. 

An important choice is thus the calibration of the regularization parameters Xj. To illustrate the 
benefits of using the smoothed Frechet mean with an automatic choice of the Aj's by GOV, we used 
M = 100 replications of model ()2.2|) with J = 15. For each repetition m = 1, . . . ,M, we have computed 
an estimator fp^^lcmstes °^ f using the Procrustes mean without smoothing, and an estimator f^l^j^^^ of 
/ using a smoothed Frechet mean where each Xj is chosen by denoising the j-th signal by GCV. In Figure 
[71 we display boxplots (for m = 1, . . . , M) of the mean squared errors 



n 

=1 £=1 



It can be seen that, over M = 100 replications, the mean squared error of the smoothed Frechet mean is 
clearly much lower than the one obtained by Procrustes mean. 
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(e) (f) 

Figure 6: (a) Unknown mean shape /, (b)-(d) Three noisy signals out of J = 15 sampled from model 
(12. 2p . (e) Estimation of / by smoothed Frechet mean (solid curve) using low-pass Fourier filtering with 
Xj chosen by GCV, (f) Estimation of / by Procrustes mean without smoothing (solid curve). The dashed 
curve in (e) and (f) is the mean shape /. 



Figure 7: Boxplots of MSE{fp^l^^J (left-hand boxplot) and M5S(/)I"j^^^„^^^J (right-hand boxplot) for 
m = 1, . . . ,M. 



p(m) 
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3 Application to ECG data analysis 



Now, let us return to the analysis of the ECG records displayed in Figure [TJ We propose to compare the 
results obtained by the smoothed Frechet mean using either translation operators or non-rigid operators 
given by the family of parametric diffeomorphisms of [0, 1] described in Section 12.41 The smoothing of 
the J signals obtained after segmentation of the ECG records is done by wavelet thresholding with a 
data-based choice of the regularization parameters dj as explained in Section [2.21 

The computation of smoothed Frechet means using translation operators in the normal and cardiac 
arrhythmia cases are displayed in Figure [8] and Figure [H The results obtained using the Procrustes mean 
are very similar, and they are thus not reported. Note that using Frechet or Procrustes means with 
translation operators is very similar to the alignment of signals using cross-correlation which is the widely 
used technique for signal averaging in ECG data analysis. 

One can see that, for the normal case, the results are very satisfactory as the mean heart cycle displayed 
in Figure [8^a) is a good estimation of the typical shape of the signals displayed in Figure [8l^b)-(g). For the 
case of cardiac arrhythmia, the results are not very satisfactory since there is still a low-pass filtering effect 
on the shape of the QRS complex (around the time point t ~ 0.45) in the mean heart cycle displayed in 
Figure m^a). In this case, the Frechet mean does not represent very well the shape of the signals in Figure 
[Hb)-(g). This low-pass filtering effect is due to a much larger variability in shape of the signals (after 
segmentation) in the case of cardiac arrhythmia than in the normal case. Since the activity of the heart 
can be very irregular in cases of cardiac arrhythmia, modeling shape variability using only translation 
operators is not flexible enough. A more precise alignment to take into account a local variability in lag 
and duration of the heart beats is thus needed. 

The computation of smoothed Frechet means using non-rigid operators in the normal and cardiac 
arrhythmia cases are displayed in Figure [10] and Figure [TTJ In the normal case, the result is very similar 
to the one already obtained by the Frechet mean using translation operators, as shown by the comparison 
of Figure [8] and Figure [TOl In the case of cardiac arrhythmia, the results displayed in Figure [9] and Figure 
[TT] clearly show the improvements obtained by using non-rigid deformation operators. The shape of the 
mean heart cycle in Figure [TlTa) is a better estimation of the typical shape of the signals displayed in 
Figure [T]T b)-(g). In particular, there is no low-pass filtering effect around the time point t ~ 0.45 in the 
shape of the Frechet mean. 

4 Discussion and conclusion 

We have presented a new algorithm for aligning heart beats extracted from an ECG record. Our ap- 
proach is based on the notion of smoothed Frechet means of curves using either translation or non-rigid 
deformation operators. In the case of a normal ECG obtained from a healthy subject, our approach 
yields a satisfactory mean heart cycle. Note that using Frechet mean with translation operators is very 
similar to the computation of a mean heart cycle by aligning the data using cross-correlation which is 
the widely used approach in ECG data analysis for signal averaging. When using non-rigid operators to 
align hearts beats having a high shape variability, with peaks showing an important variability in lag and 
duration from one pulse to another, our approach yields significant improvements over signal averaging 
by cross-correlation. The benefits of our procedure has been demonstrated for an ECG recording of a 
subject showing evidence of significant arrhythmia. We hope that the methods presented in this paper 
will stimulate further investigation into the development of better alignment procedures that take into 
account local shape variability in heart beats extracted from ECG records. 
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(a) (b) (c) (d) 




(e) (f) (g) 



Figure 8: Normal case: (a) Frechet mean using translation operators of the J = 93 signals after segmen- 
tation of the ECG record, (b)-(g) Superposition of six signals containing a single QRS complex (dashed 
curves) with the Frechet mean using translation operators (solid curve). 




(e) (f) (g) 



Figure 9: Case of cardiac arrhythmia: (a) Frechet mean using translation operators of the J = 72 signals 
after segmentation of the ECG record, (b)-(g) Superposition of six signals containing a single QRS 
complex (dashed curves) with the Frechet mean using translation operators (solid curve). 
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(a) (b) (c) (d) 




(e) (f) (g) 



Figure 10: Normal case: (a) Frechet mean using non-rigid operators of the J = 93 signals after segmen- 
tation of the ECG record, (b)-(g) Superposition of six signals containing a single QRS complex (dashed 
curves) with the Frechet mean using non-rigid operators (solid curve). 




(e) (f) (g) 



Figure 11: Case of cardiac arrhythmia: (a) Frechet mean using non-rigid operators operators of the J = 72 
signals after segmentation of the ECG record, (b)-(g) Superposition of six signals containing a single QRS 
complex (dashed curves) with the Frechet mean using non-rigid operators (solid curve). 
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